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RASTERIZATION OF LINES IN A CYLINDRICAL VOXEL GRID 

This invention relates to computer graphics, and in particular, to the rendering 
of a straight line on a display in which the voxels are arranged on a cylindrical 
coordinate grid. 

BACKGROUND 

Many computer-implemented displays consist of two-dimensional arrays of 
individual picture elements, or pixels. To form an image, a rasterizer selectively colors 
the pixels. Because the individual pixels are so small, the display appears to a human 
viewer to be a continuous rendering of an image. This illusion is particularly effective 
for complex images of continuous tones such as photographs. 

For simple geometric shapes, however, the pixelated nature of the display can 
become apparent to the human viewer. For example, if the rasterizer is instructed to 
draw a straight line, there is no guarantee that the points on that desired line will 
coincide with the pixels that are available for rendering the line. As a result, the desired 
line is often rendered by a rasterized line of pixels that are close to but not necessarily 
coincident with the desired line. This results in rasterized lines that have a jagged or 
echeloned appearance. 

In the course of rendering an image, a large number of straight lines and line 
segments are often drawn. As a result, given a desired line, the rasterizer must 
frequently select those pixels that will minimize the jagged appearance of the resulting 
rasterized line, A straightforward mathematical approach is to use the equation of the 
desired line and the coordinates of the available pixels to minimize a least square error 
across all points on the line. While such an approach has the advantage of globally 
optimizing the selection of pixels on the rasterized line, the large number of floating 
point operations required causes this approach to be prohibitively time-consuming. 

To meet constraints on speed, rasterizers typically implement rasterization 
methods that avoid time-consuming floating point operations. Among the methods that 
meet the foregoing constraints is that taught in Bresenham> "Algorithm for 
Computer Control of a Digital Plotter," IBM System Journal, Vol. 4, pp.25-30, 1965, 
the contents of which are herein incorporated by reference. Using only integer 
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5 operations, the Bresenham algorithm reduces the choice of what pixel to select to an 
examination of the sign of a discriminant. 

Even faster rasterization methods exist that select multiple points on the 
rasterized line based on the outcome of a single decision. There also exist a variety of 
rasterization methods aimed at rendering conic sections and quadric sections on a 
1 0 pixelated display. 

However, the foregoing rasterization methods all rely on the assumption that the 
array of pixels is arranged in a uniform rectangular grid that can readily be modeled by 
a Cartesian coordinate system. This is a reasonable assumption given the prevalence of 
two-dimensional displays such as computer monitors and printers at the time these 
1 5 algorithms were developed. 

Since then, however, volumetric, or three-dimensional displays have been 
developed. Such displays permit the generation, absorption, or scattering of visible 
radiation from a set of localized and specified regions within a volume. Examples of 
such systems are taught in Hirsch U.S. Patent No. 2,967,905, Ketchpel U.S. Patent No. 
20 3,140,41 5, Tsao U.S. Patent No. 5,754,147 and on pages 66-67 of Aviation Week, 
October 31, 1960. 

In a typical volumetric display 1, shown in FIG. 1 , a motor 2 spins an imaging 
plate 3 rapidly about an axis 4. A light source 5 under the control of a rasterizer 6 
illuminates selected spots 7 on the imaging plate 3 at successive instants. If the imaging 
25 plate 3 spins rapidly enough, and if the successive instants are separated by short 
enough time intervals, a continuous curve will appear to hang in mid-air. 

FIG. 2 illustrates the projection 8 onto a plane perpendicular to the axis 4 of an 
exemplary line rendered by the rasterizer 6. In FIG. 2, the imaging plate 3 is shown in 
six successive instants as it rotates around the axis 4, now perpendicular to the page. At 
30 each of the six instants, the light source 5, under the control of the rasterizer 6, 

illuminates a spot 7 on the imaging plate 3. As shown in FIG. 2, by aiming the light 
source 5 at the correct spot and firing the light source 5 at the right time, it is possible to 
trace out the line. It is the function of the rasterizer 6 to correctly aim and fire the light 
source 5 so as to trace out the line. 
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To aim the light source 5, the rasterizer 6 needs a way to identify points in 
space. In other words, the rasterizer 6 needs a coordinate system. One possible 
coordinate system is a Cartesian coordinate system. 

Using a Cartesian coordinate system, the rasterizer 6 would draw the line 8 in 
FIG. 2 by specifying, for example, that at time t\ 9 the light-source 5 should aim 30 units 
north, at time t 2 , the light-source 5 should aim 29 units north, 1 unit east and 1 unit in 
elevation, at time the light-source 5 should aim 28 units north, 2 units east and 2 
units in elevation, and so on. To ensure that the imaging plate 3 is at the appropriate 
location, the rasterizer 6 also has to specify the values of the times t\ 9 1 2 and h based on 
how fast the imaging plate 3 is spinning. Although the spinning of the imaging plate 3 
can be resolved into a north-south component and an east-west component, this is a 
computationally taxing exercise that can easily be avoided by using a polar coordinate 
system. 

As an alternative, the rasterizer 6 could draw the line in FIG. 2 in cylindrical 
coordinates by specifying, for example, that at time t\ 9 the light source 5 should aim 30 
units away from the axis 4 at a 90 degree angle, at time / 2 , the light source 5 should aim 
29.02 units from the axis 4 at an angle of 88.03 degrees and point up by 1 unit of 
elevation, and at time / 3 , the light source 5 should aim at 28.07 units from the axis 4 at 
an angle of 85.91 degrees and point up by 2 units of elevation. This would, of course, 
draw the same line that was drawn using Cartesian coordinates. The advantage of using 
cylindricaT coordinates becomes apparent when it comes time to specify when to fire 
the light source 5. Because the imaging plate 3 is spinning, it is far more natural to 
represent its motion in terms of degrees per second than it is to resolve its motion into a 
north-south component and an east-west component. For this reason, calculations 
involving motion of the imaging plate 3 in a volumetric display are best performed in 
terms of a cylindrical coordinate system. 

The process of rendering a line in a cylindrical coordinate system is 
fundamentally different from that of rendering a line in a Cartesian coordinate system. 
In a Cartesian coordinate system, a line has a constant slope. A unit change in one 
direction always results in the same change in the remaining two directions, regardless 
of where that change occurs. This is not the case in a cylindrical coordinate system. 
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5 In a cylindrical coordinate system, the slope of a line can vary dramatically with 

position along the line. In contrast to a Cartesian coordinate system in which the slope 
of a line is constant at all points on the line, the slope of a line in a cylindrical 
coordinate system can vary significantly along the line. This is because the coordinate 
grid in a cylindrical coordinate system is not a spatially uniform grid, as is the 
10 coordinate grid in a Cartesian coordinate system. 

For those portions of the line that are closest to the axis of the cylindrical 
coordinate system, the change in distance of the line from the origin changes very 
slowly with angle. For those portions of the line that are far from the axis of the 
cylindrical coordinate system, even a small change in angle results in a huge change in 
1 5 radial distance from the axis. 

SUMMARY 

A rasterizer according to the invention recognizes that the dominant direction of 
a desired line in a cylindrical coordinate system changes in a predictable manner. The 
rasterizer uses this predictable sequence of dominant directions, referred to as a 

20 "dominance progression," to avoid having to determine the dominant direction at each 
point on the desired line before rendering a voxel representative of that point Instead, 
the rasterizer detects boundaries between constituent line-segments of the desired line 
and uses the dominance progression to determine the dominant direction of the line at 
any point Because the process of detecting a transition between boundaries is much 

25 faster than that of evaluating a dominant direction, a rasterizer practicing the method of 
the invention can rapidly and efficiently draw a rasterized line in a cylindrical 
coordinate system. 

The rasterizer renders, on a volumetric display, a rasterized line that 
approximates a desired line. It does so by rendering, on the volumetric display, a first 
30 voxel representative of a first point from the desired line. The first point is selected 
from a first constituent line-segment having a first dominant direction. Then, the 
rasterizer selects a second point from the desired line and determines whether that 
second point is on a second constituent line-segment, with a second dominant direction, 
or whether that second point is on the same constituent line segment. 
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5 On the basis of this determination, the rasterizer renders a second voxel 

representative of the second point on the desired line. This second voxel is located so as 
to minimize a distance between the second voxel and the second point on the desired 
line. The second voxel can be selected, for example, by application of a midpoint 
method to select the second voxel from a plurality of candidate voxels. 

10 In one practice of the invention, the rasterizer evaluates a difference between a 

threshold and a transition-indicating expression evaluated at the second point. This 
difference is indicative of a transition from the first constituent line-segment to the 
second constituent line-segment The transition indicating-expression is typically an 
expression whose value depends on a slope of the desired line at the second point. The 

1 5 slope can be defined as a ratio of how much the desired line changes in the first 

dominant direction with respect to a change in one of two non-dominant directions. 

In some cases, the orientation of the line to be rasterized is such that the 
transition cannot be reliably detected. This occurs, in part, because of computation on a 
digital computer cannot be carried out with infinite precision. In one practice of the 
20 invention, erroneous transitions are avoided by selecting the transition-indicating 

expression to be one that identifies a line having a point at which the slope is equal to 
the threshold. 

Unless otherwise defined, all technical and scientific terms used herein have the 
same meaning as commonly understood by one of ordinary skill in the art to which this 

25 invention belongs. Although methods and materials similar or equivalent to those 

described herein can be used in the practice or testing of the present invention, suitable 
methods and materials are described below. All publications, patent applications, 
patents, and other references mentioned herein are incorporated by reference in their 
entirety. In case of conflict, the present specification, including definitions, will control. 

30 In addition, the materials, methods, and examples are illustrative only and not intended 
to be limiting. 

Other features and advantages of the invention will be apparent from the 
following detailed description, and from the claims. 



5 



WO 02/33661 



PCT/USO 1/32246 



5 BRIEF DESCRIPTION OF THE DRAWINGS 

FIG. 1 shows a volumetric display for practice of the rasterization method of the 
invention; 

FIG. 2 is a projection of a desired line on the r-9 plane; 

FIGS. 3-4 show two views of the desired line and a cylindrical surface to which it 
10 is tangent; 

FIGS. 5, 6, 7A and 7B show curves derived from three discriminant functions used 
in rasterizing the desired line; 

FIG. 8 is a flow chart summarizing the steps in the practice of the invention; 

FIGS. 9-1 1 show curves derived from the three discriminant functions of FIGS. 6, 
15 6, and 7A together with the locus of points for which their derivatives are equal 

to one; 

FIGS. 12-13 are examples of cases in which the proximity of a discriminant 
function to the locus of points in FIGS. 10 and 1 1 gives rise to numerical 
instability; 

20 DETAILED DESCRIPTION 

A line traversing a cylindrical coordinate system defines a sequence of points, 
each of which is defined by a radial coordinate r, an angular-coordinate 0, and an axial 
coordinate h. As one traverses the line, these points change by different amounts 
depending on the direction of the line. For example, if the line is nearly parallel to the 

25 axis of the cylindrical coordinate system (hereafter referred to as the "axis 5 '), r and 0 
will change by only small amounts compared to the change in Such a line is referred 
to as "^-dominant" because its dominant direction is the axial direction. If the line 
passes through the axis of the cylinder and deviates only slightly from a plane 
orthogonal to that axis, then 0 will never change at all and h will change by only a 

30 small amount. Such a line is referred to as "r-dominant" because its dominant direction 
is the radial direction. For short lines that do not pass through the axis, 0may change 
more rapidly than either r or h. Such lines are said to be "^dominant" because their 
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5 dominant direction is the angular direction. 

As a general rule, the dominant direction of a line in a cylindrical coordinate 
system changes as one traverses the line. It is useful therefore to consider a line as a 
concatenation of line-segments, each of which has a dominant direction that differs 
from its neighboring line-segments. Within each line-segment, the dominant direction 
10 of the line is the same. However, when crossing the boundary between two line- 
segments, the dominant direction changes. 

The geometric meaning of "line" as a set of points extending infinitely in both 
directions is not useful in a discussion of a volumetric display having finite extent. 
Hence, throughout this specification, the word "line" generally refers to line-segments 
1 5 defined by two endpoints. With this definition, it becomes meaningful to refer to lines 
as being "short" or "long." A "short line" is a line-segment whose length is short 
relative to the extent of the volumetric display. A "long line" is a line-segment other 
than a short line. 

In the rasterization method of the invention, the rasterizer 6 begins at one 
20 endpoint of the line and identifies the dominant direction of the line-segment containing 
that endpoint To reach the next point on the line, the rasterizer 6 first increments the 
value of the coordinate (hereafter referred to as the "dominant coordinate) associated 
with the dominant direction for that line-segment. Then, the rasterizer 6 decides 
whether to increment the values of coordinates associated with one or the other non- 
25 dominant direction. 

To determine whether to increment the values of the coordinates associated with 
the non-dominant directions, it is useful to define three discri m ina n t functions, each of 
which defines a surface in the cylindrical coordinate system. Each of the three 
coordinates is an argument to two of these discriminant functions. The discriminant 
30 functions are selected so that the intersection of any two of them defines the desired 
line. For convenience of notation, these discriminant functions are referred to as: 

F(r,0) = O 

G(r, 70 = 0 
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The functional forms of the above three discriminant functions are best 
understood with reference to FIGS. 3 and 4. As shown in those figures, a desired line 
10 in cylindrical coordinates can be conveniently described in terms of four geometric 
measurements. A tangent plane containing the desired line 10 is defined by the line's 
10 point of tangency 14 (r n a,h t ) with a cylindrical surface 15. The specific line within 

this tangent plane is defined by the angle of twist, /?, of the desired line 10 in the 
tangent plane. 

The first discriminant function, F(r 9 0) , is the projection of the desired line 10 
onto a plane that is perpendicular to the axis of the cylindrical coordinate system: 

F(r,0) = reos(0-a)-r, =0 

where r t is the radial distance to the point of tangency 14 and a is the angular 
coordinate at the point of tangency 14. An intersection of the r-6 plane and the surface 
defined by F(r 9 0) = 0 is shown in FIG. 5. 

The second discriminant function, G(r 9 h) , is a positive branch of a hyperbola. 
The hyperbola has an axis of symmetry extending in the h direction and passing 
through the point of tangency 14 at (r h hi). The angle formed by the r-axis and the 
hyperbola's asymptotes is /?. The hyperbola is derived from the desired line 10 using 
the equation: 

G(r 9 h) = Tj(r 2 -r t 2 )-(h-h t ) 2 =Q 

where rj is defined as tan 2 ((3). FIG. 6 shows the curve satisfying the constraint 
G{r 9 h) = 0 in the r-h parametric domain. 

The third discriminant function, H(h 9 0) , defines a tangent function in the 9-h 
parametric domain. The inflection point of this tangent function is located at the point 
of tangency 14, (0 = a 9 h = h t ) . The equation is given by: 
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5 H(h, 0) = f tan(0 -a) + h t -h = Q 

where £is defined to be r t tan(yff) . FIGS. 7 A and 7B show curves defined by 
H(h 9 0) = 0 in the O-h parametric domain corresponding to positive and negative 
values of £ 

Note that unlike the first two discriminant functions, the sign of the slope of 
1 0 H(h,0) changes depending on the sign of £ In implementing the third discriminant 
function, the rasterizer 6 can use two separate (but very similar functions). 
Alternatively, the rasterizer 6 can incorporate the sign of £ into a single function call. 
The use of two separate functions may be more efficient because doing so avoids extra 
operations associated with changing signs of subexpressions. 

15 A dominant direction characterizes each point on the desired line 10. The 

dominant direction associated with a point on the desired line is that direction in which 
the coordinate for that direction changes most rapidly relative to the other coordinates. 
The rasterizer 6 identifies this direction by using the above three functions to evaluate 
the derivatives dr/dh , dr/dO , and dh/d9 . Because of the constraints 

20 F(r 9 0) = G(r,h) = H(h,0) = 0 , it follows that dF = dG = dH = 0. Hence, one obtains 
the following expressions: 

dr , n .radians 

rtan(0 — a) — 

d0 step 

dr _ h-h t 
dh rjr 

dh K 2 , n .radians 

— = £ sec (0 - a) 

d0 step 

25 To improve performance, angular coordinates are expressed as dimensionless 

integers rather than in radians. A "radians/step" constant for converting between 
dimensionless integers and radians is therefore incorporated within the 0 argument of 
the trigonometric functions. The scale factor "radians/step" that appears in the first and 
third derivatives thus arises from the chain rule of elementary calculus. 
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5 Given a point on the desired line 10, one can evaluate the three discriminant 

functions F(r,0), G(r,h) and H(h,0) . The derivatives of these three functions at that 
point on the desired line 10 identify the dominant direction at that point. Once the 
dominant direction is known, there exist known midpoint methods for determining 
successive voxels, thereby generating a rasterized line that approximates the desired 
10 line 10. 

Because of the time required to evaluate the three discriminant functions and 
their derivatives at each point on the desired line 10, it is preferable to instead predict 
the progression of dominance one encounters while progressing along the desired line 
10 and to identify indicia of a change in the dominant direction of the desired line 10. 
15 For a line in cylindrical coordinates, a dominance progression can be derived by 
observing the behavior of the three discriminant functions as 0 proceeds from 
0 = a-7v/2 to 6 = a + 7t/2 . 

The dominance progression for the discriminant function F(r,0) is that of a 
secant function. As 0 increases, F(r,0) , and hence the desired line 10 from which it 
20 was derived, progresses from being r-dominant, to being ^-dominant, and then back to 
being r-dominant. 

The dominance progression for the discriminant function G(r,h) is that of one 
branch of a hyperbola. In this case, the dominance progression depends on the slopes of 
the asymptotes. If the slopes of the asymptotes have a magnitude greater than one, then 

25 G(r 9 h) , and hence the desired line 10 from which it was derived, progresses from 

being r-dominant, to being ^-dominant, and then back to being r-dominant, switching 
sign at the point of tangency 14 (h t ,r i ) . However, if the slopes of the asymptotes have 
magnitudes that are less than one, the hyperbola never becomes r-dominant. Under 
these circumstances, h is more dominant than r at all points on the desired line 10. 

30 Geometrically this corresponds to a desired line 10 that is nearly parallel to the 
cylindrical axis. 

The third discriminant function, the tangent function H(h 9 0) , generally 
progresses from being /z-dominant, to being ^-dominant and finally back to being h- 

10 
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5 dominant. However, the minimum slope of a tangent function is equal to its coefficient, 
£ Hence, if <£is greater than one, then the slope of the tangent function H(h,0) is 
everywhere greater than one. Thus it is possible for a line to be A-dominant everywhere 
and never ^dominant Geometrically this corresponds to a desired line 10 for which 
either r t is large, or to a desired line 10 that is nearly parallel to the cylindrical axis. 

10 The dominance progression for each of the three discriminant functions can 

therefore be summarized in the following table: 



increasing 
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r 


G(r,h) 


Or) 


h 


to 


H(h,6) 


h 


(0) 


h 



In the above table, parentheses enclose those coordinates that, for certain lines, may 
never become dominant. Note that the inflection point on each of these three 
15 chscriminant functions coincides with the line's point of tangency 14 (r, ,a, h t ) with the 
cylindrical surface 15. 

The middle column of the table indicates that at the point of tangency 14, r can 
never be dominant. At that point, only 0 or h can be dominant. Which of the two is 
dominant can be determined by observing the behavior of H(h 9 0) . The dominance 
20 progression for H(h,0) indicates that if there exists a segment on the desired line 10 in 
which 0 is more dominant than h 9 the point on that line at which 0 = a will be in that 
region. Consequently, 0 dominance, if present at all, will manifest itself near the 
6 = a point on the line. 

Proceeding outward from the 0 = a point on the desired line 10, the next 
25 coordinate to become dominant will be either r or h. The dominance progression of the 
hyperbola, G(r,/z) indicates which one it will be. Based on the dominance progression 
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of G(r 9 h) , it is apparent that if there exists an /z-dominant line-segment on the desired 
line 10, that line-segment will be closer to the point of tangency 14 than any r-dominant 
segment. 

Based on the foregoing, it is always the case that, in traversing a line in 
cylindrical coordinates from 6 = a-n/2 to & = a + x/2 , the dominant directions of 
the consecutive line-segments that constitute the desired line 10 will progress in the 
following order: 

Note that the above dominance progression does not require that each direction 
be dominant at some point on the desired line 10. Instead, the above dominance 
progression states that if the desired line 10 includes a line-segment in which a 
particular direction dominates, then that line-segment will appear in relation to the 
other line-segments as indicated above. In referring to a line-segment, a or 6 — " sign 
indicates whether that line-segment lies before the point of tangency 14 or after the 
point of tangency 14. Thus the constituent line-segments of the desired line 10 are 
conveniently labeled as follows: 

The recognition of the foregoing dominance progression in any line in 
cylindrical coordinates enables the implementation of a rasterization method that does 
not require computationally expensive evaluation of the three discriminant functions to 
identify the dominant direction at each point on the line. Instead, rasterization proceeds 
following steps shown in FIG. 8. 

The rasterization method shown in FIG. 8 begins with an initialization step in 
which the rasterizer 6 renders the first voxel (step 16). This first voxel corresponds to a 
first endpoint of the desired line 10. The first endpoint is selected so that traversing the 
desired line 10 from the first endpoint to the second endpoint results in a decreasing r 
coordinate followed by an increasing r coordinate. 

The initialization step is followed by steps in which the rasterizer 6 proceeds 
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from one constituent line-segment of the desired line 10 to the next. Within a 
constituent line-segment, the rasterizer 6 steps unconditionally in the dominant 
direction of the line-segment. It then determines whether or not to also step in the non- 
dominant directions. To do so, the rasterizer 6 evaluates a discriminant function (one of 
F(r,0) , G(r,h) , or H(h 9 0) ) associated with each of the other two non-dominant 
directions. After rendering each voxel, the rasterizer 6 updates all the slope values for 
the discriminant functions and checks to see if the next point is within the same line- 
segment. 

Following the initialization step, the rasterizer 6 renders voxels representative of 
points in the R~ line-segment of the desired line 10, if any (step 18). In doing so, the 
rasterizer 6 evaluates the signs of the first and second discriminant functions, F(r,0) 
and G(r,h) , to determine whether to step in the 0and h directions respectively. Before 
rendering a voxel representative of a point on the desired line 10, the rasterizer 6 
verifies that that point is still in the RT line-segment of the desired line 10. If it is not, 
then the rasterizer 6 moves to the step of rendering voxels representative of points on 
the H ' line-segment of the desired line 10, if any (step 20). 

Once the rasterizer 6 determines that it has crossed over into the H " line- 
segment of the line, it evaluates the signs of the second and third discriminant 
functions, G(r 9 h) and H(h,0) , to determine whether to step in the r and 0 directions 
respectively. Before rendering a voxel representative of a point on the desired line 10, 
the rasterizer 6 checks to see if that point is still in the H~ line-segment of the desired 
line 10. If it is not, then the rasterizer 6 moves to the step of rendering voxels ~, 
representative of points on the 0~ line-segment of the desired line 10, if any (step 22)7 

Once the rasterizer 6 determines that it has crossed over into the 9-line-segment 
of the line, it evaluates the signs of the first and third discriminant functions, F(r 9 0) 
and H(h,0) , to determine where to step in the r and h directions respectively. The 
rasterizer 6 continues to operate in this manner as long as 0 remains the dominant 
direction. Hence, the rasterizer 6 operates in this manner even after crossing over into 
the 0* line-segment of the desired line 10 (step 24). 
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Once in the 0* line-segment of the desired line 10, the rasterizer 6 checks to see 
if the next point to be drawn is still in the# + line-segment of the desired line 10. If it is 
not, then the rasterizer 6 moves to the steps of rendering voxels representing points on 
the H + line-segment of the line (step 26), if any, and voxels representing points on the 
R + line-segment of the line (step 28), if any. These steps are carried out in a manner 
analogous to the earlier steps of rendering voxels representing points in the H~ and R~ 
line-segments of the desired line 10. 

To locate the boundaries of the constituent line-segments of the desired line 10, 
the rasterizer 6 determines the slope of the discriminant function with respect to the 
currently dominant variable. If this slope transitions from being less than one to being 
greater than one (or vice versa), the rasterizer 6 interprets this as indicating that it has 
now crossed into a different constituent line-segment of the desired line 10 and that the 
dominant direction has changed in accordance with the progression of dominance 
identified above. 

FIG. 9 shows the curve F(r 9 0) = 0 (the solid curve) together with the locus of 
all points for which dr/dO is equal to one (the dashed curve). As is apparent from the 
figure, there are only two points in r-8 space that are both on the discriminant function 
and for which drjdO is equal to one. These two points, which correspond to the 
intersections of the dashed curve and the solid curve in the figure, unambiguously 
identify the points at which the line transitions from one constituent line-segment to the 
next. FIGS. 10 and 1 1 show the curves G(r,/z) = 0 and H(h 9 0) = 0 (solid curves) and 
the corresponding loci of points for which dr/dhand dhjdO are equal to one. As 
shown in FIGS. 10 and 1 1, there are only two points in each case that are both on the 
discriminant function and for which the derivative is equal to one. 

For each of the three cases above, there also exist points that almost, but do not 
quite satisfy the two conditions set forth for each case. These points correspond to those 
that are in the neighborhood of the intersection of the dashed and solid curves in FIGS. 
9-1 1 . One such set of erroneous points is illustrated by the shaded region in FIG. 10. 

Because the only points that can be rendered on the display are those that fall on 
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the cylindrical coordinate grid, it is possible that the rasterizer 6 will mistakenly 
consider one of the points in the set of erroneous points to be the actual intersection of 
the dashed and solid curves. This may result in the rasterizer 6 assuming that a 
transition from one constituent line-segment to the next has occurred when in fact, no 
such transition has occurred. 

For a desired line 10 whose discriminant functions are as shown in FIGS. 9-11, 
a numerical error of the type described above will have little or no practical 
consequences. This is because in each of FIGS. 9-1 1, all points in the set of erroneous 
points are proximate to the true intersection point of the dashed and solid curves. 

For some lines in a cylindrical coordinate system, however, the set of erroneous 
points can extend further away from the true intersection point of the dashed and solid 
curves. For example, referring back to FIG. 4, when the angle of twist p in the tangent 
plane approaches 45 degrees, tan(j3) approaches one. When this occurs, the 
discriminant function G(r 9 h) and the locus of points for which drjdh is equal to one 
become close to each other even at points that are far away from their intersection. In 
this case, the set of erroneous points, shown in FIG. 12, extends far beyond the 
intersection of the two curves. Under these circumstances, an erroneous selection of 
one of these points can results in a grossly premature transition from one constituent 
line-segment to the next. 

Another example of a desired line TO that may result in a numerical error by the 
rasterizer 6 is one that passes very close to the cylindrical axis (r, « 0) or one in which 
the desired line 10 lies mostly in the r-theta plane (J3 « 0) . In either of these cases, 
£ approaches zero. This causes the vertical dashed lines of FIG. 1 1 to move outward 
toward the tail of the tangent function, as shown in FIG. 13. Under these circumstances, 
the vertically extending set of erroneous points can include erroneous points are far 
away from the actual intersection of the dashed lines and the solid curve. This again can 
result in an inappropriate transition from one constituent line-segment to the next. 

One way to eliminate numerical errors of the type described above is to 
symbolically compute the coordinates of the actual intersection of the dashed and solid 
curves in FIGS. 9-11 and to use those coordinates as transition thresholds. In practice it 
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5 suffices to compute only one of the two coordinates that define the intersection. Of the 
two coordinates that define the intersection point, the more reliable is the dominant 
coordinate. For example, when in an r-dominant line-segment the rasterizer 6 would 
compare the current value of r on the desired line 10 with the r coordinate of the 
intersection to determine whether to transition out of the r-dominant line-segment. In 
1 0 terms of FIGS. 9-1 1 , the symbolically computed coordinate defines either a horizontal 
or vertical line which goes through the intersection point. 

Using the equations for F(r 9 0) and dr I d6 , one can eliminate a common 
variable r and derive the following expression for a vertical line in FIG. 9: 

radians 

Fi(0) = cos 2 {G-d)- r t sin(0 - d) 

step 

15 Similarly, using the equations for G{r y h) and dr/dh, one can eliminate the common 
variable 0and derive the following equation for a horizontal line in FIG. 10: 

G/(r) = (l-^)r 2 -r f 2 
Finally, using the equations for H(h,0) 3ixddH/d0 , one can again eliminate the 
common variable 0 and derive the following expression for a horizontal line in FIG. 1 1 : 

20 Hm = (h-h t ) 2 - ihSgn - 

rad 

These three transition-indicating expressions are used as substitutes for 
expressions that are prone to numerical instability under the circumstances set forth 
above. For example, if the rasterizer 6 is in the R~ line-segment of a desired line 10, it 
can recognize a transition to a successor line-segment upon recognition of the truth of 
25 either one of the following Boolean expressions: 

dr t 

— > — 1 , or 

d0 

— • hSgn > -1 
dh 

However, the second expression is prone to numerical instability for certain lines. As a 
result, the Boolean expression 

30 ~hSgn>-l 
dh 
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is replaced by a Boolean expression that uses the more stable transition-indicating 
expression Gi(r) : 



10 



G/(r)<0. 

The following table lists, for each constituent line-segment, a pair of Boolean 
expressions that can be used to indicate a transition from one constituent line-segment 
to another: 



Constituent line segment 


Expression 1 


Expression 2 


R~ 


dr 

— >-l 
d& 


^-•hSgn > -1 
dh 


H~ 


dh , 0 
~hSgn<\ 


dr , _, 

— hSgn<0 

dh 


0- 


le >0 


none 


e + 




dr 

— >1 
d9 




~hSgn>\ 
dh 


none 




none 


none 
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Certain of the foregoing Boolean expressions are prone to numerical instabilities. To 
avoid numerical instabilities, these Boolean expressions can be replaced by more stable 
Boolean expressions involving the transition-indicating expressions Fi, Gi and Hi. The 
following table lists alternative pairs of Boolean expressions for determining the 
occurrence of a transition between constituent line-segments. In the following table, all 
Boolean expressions that are prone to numerical instability have been replaced by more 
stable expressions involving the transition-indicating expressions Fi 9 Gi and Hi. 



Constituent line segment 



Expression 1 



Expression 2 
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dr 

Te > ~ 1 


Gi<0 


H~ 


Hi<0 


h-hSgn > h t hSgn 




da 


none 




~hSgn>\ 
uu 


Fi<0 




Gi>0 


none 


R* 


none 


none 



5 

The code provided as an appendix below shows one implementation of a 
rasterizer according to the invention. However, it is to be understood that while the 
invention has been described in conjunction with the detailed description thereof, the 
foregoing description and the accompanying code is intended to illustrate and not limit 
10 the scope of the invention. 



APPENDIX 



void writeLineFloatCyIindrical(float xl, float yl, float zl, 
float x2, float y2, float z2) 

15 // THE FOLLOWING STATIC CONST VALUES ARE FIXED FOR A GIVEN DEVICE 

static const float radsPerStep = M_2PI/g_stepsA; 

static const float stepsPerRad = 1.0/radsPerStep; 

static const float thetaStep = radsPerStep; 

static const float halfThetaStep = 0.5*radsPerStep; 
20 static const float thetaDomLimit = 

2*IntegerDoubleSin( 1 )/(5-4*IntegerDoubIeCos( 1 )); 

static const float MACHINE JEPSILON = le-7; 
#ifdef INCREMENTAL_TRIG 

// COMPUTE VALUES FOR INCREMENTAL SIN/COS EVALUATION 
25 // aW,b W for whole steps; aH,bH for half steps 

static const float sinQuarterThetaStep = sin(halfThetaStep*0-5); 

static const float aH = 2*sinQuarterThetaStep*sinQuarterThetaStep; 

static const float bH = IntegerDoubleSin(l); 

static const float aW = 2*bH*bH; 
3 0 static const float bW = IntegerDoubleSin(2); 
^endif 

#ifdef INCREMENT AL_TRIG 
# define TRIG_HALF_STEP(c,s) \ 
35 do { prevCos = c; prevSin = s; \ 

c = prevCos - (aH*prevCos + bH*prevSin);\ 



18 



WO 02/33661 



s = prevSin - (aH*prevSin - bH*prevCos); } while(O) 

# define TRIG_WHOLE_STEP(c,s) \ 
do { prevCos = c; prevSin = s; \ 

c = prevCos - (aW*prevCos + bW*prevSin);\ 

s = prevSin - (aW*prevSin - bW*prevCos); } while(0) 

#else 

# defi ne TRIG_HALF_STEP(c,s) \ 

do { c - cos( m_th*radsPerStep - tht ); \ 

s = sin( m_th*radsPerStep - tht ); } while(0) 

# define TRIG_WHOLE_STEP(c,s) \ 

do { c - cos( m_th*radsPerStep - tht ); \ 

s = sin( m_th*radsPerStep - tht >; > whiJe(O) 
#endif 

///////////////////////////// LINE SETUP IllllllltlllllllllllWtllimil 

II Test ordering and ensure counterclockwise orientation 
bool ccw = (xl *y2-x2*yl) > 0; 
if (!ccw) { 

// SWAP ENDPOINTS! 
float tmp; 

tmp m xl ; xl=x2; x2 = tmp; 
tmp = yU yl^; y2 = tmp; 
tmp = zl; z\ s=s z2\ z2 = tmp; 

} 

float dx = x2-xl; 
float dy^yl-yl; 
float dz = z2-zl; 
float rlf=sqrt(xl*xl + yl*yl); 
float r2f = sqrt(x2*x2 + y2*y2); 
float rl - halfround( rlf ); 
float r2 = half_round( r2f ); 
float thlf, thl,th2; 

bool thlgnore = false; // Default is that theta could be dominant 
// Handle degerate endpoint-on-origin case by virtually 
// perturbing origin endpoint towards other endpoint 
if(xl=0&&yl=0){ 

thlf = atan2(y2, x2) * stepsPerRad; 

thl » th2 - half_round( thlf); 

thlgnore = true; 

} 

elseif(x2=0&&y2=0){ . . . 

thlf ■ atan2(yl, xl) * stepsPerRad; 
thl =half_round( thlf); 
th2 = thl; 
thlgnore = true; 

} 

else { 

thlf - atan2(yl, xl) * stepsPerRad; 
thl =half_round( thlf); 

th2 = half_round( atan2(y2, x2) * stepsPerRad ); 
if (thl=th2) thlgnore = true; 

} 

float hi = half_round( zl ); 
float h2 = half_round( z2 ); 

// Check to see if we're going around the wrong way. Fix if so. 
if(thl>th2) 
th2 += g_stepsA; 

// COMPUTE CHARACTERISTIC LINE CONSTANTS 

// implicit fh F(r,th) = r(Acos(th)+Bsin(th)) = 0 

//const float A = dy; 

//const float B » -dx; 

//const float C = -(A*xl+B*y 1); 
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5 const int hSgn = dz<Q r >-\ : 1; 

const float hHalflnc « 0.5 * hSgn; 

float polarL2 - dx*dx + dy*dy; // length squared in 2D polar plane 
float polarL2inv, polarLinv; 
10 if (polarL2 < MACHINE_EPSILON) { 

polarL2 = 0; 

polarL2inv = 0; 

polarLinv - 0; 

15 Le{ 

polarL2inv= 1.0/poiarL2; 

polarLinv = sqrt(poIarL2inv); // inv. len. in 2D polar plane 

} 

20 const float rt = // radius at point of tangency 

(polarL2=0) ? rlf :(xl*dy-yl*dx)*polarLinv; 
const float rt2 = rt*rt; 

float ht; // height at pt of tangency 

// Use endpoint closer to polar origin for less numerical error 
25 if(rl<r2) 

ht = zl - dz*(xl *dx+yl *dy)*polarL2inv; 
else 

ht - z2 - dz*(x2*dx+y2*dy)*polarL2inv; 
const float tht = atan2(-dx,dy); // angle at point of tangency 
3 0 const float eta = (polarL2in v > 0)? dz*dz*polarL2inv : 2; 

// eta ~ 2 just convenient indication mat r never dominates h 
const float one_m_eta = 1-eta; 
const float xi = rt * dz*po!arLinv; 

35 // must decide the min radius for regions 2- and 1- 
float rrnin; 

bool start_minus;// does first endpoint lie decreasing r part of line? 
if (dx*xl + dy*yl < 0 &&. dx*x2 + dy*y2 > 0) { 

rmin = half_round(rt); 
40 if (rmin<0) rmin = -rmin; 

start_minus = (rlf >= rt); 

> 

else { 

// minimal radius appears outside the line segment 
45 rmin 88 rl<r2 ? rl : r2; 

start_minus = (rlf >= r2f); 

} 

// The maximal value for h (premultiplied by the sign of dh) 

float hmax = (hSgn*h2<hSgn*ht) ? half round(hSgn*h2) : half round(hSgn*ht); 

50 

// Because eta_expr is <0 whenever the potential region of h-dominance 
// is less than one voxel's h-width, h is never dominant, 
const float eta_expr = eta*eta*rt + eta - 1; 
^ const hlgnore = eta_expr < 0 && polarL2 != 0; 

II Theta can never be dominant outside of r=0.5 pixels 
// if rt < sin(thetastep/2); 
if (Uhlgnore) { 
thlgnore - (rt< thetaDomLimit); 

60 } 

// Left and right hand sides of slope inequality expressions 
float dhdthj, dhdthjr, 
float drdhj; 
65 float drdthj, drdthjr; 

// ... and the sign of drdth_r which is multiplied through, 
float drdth_s; 

float m_r, m_th, m_h, m_cos, m_sin; // Midpoint test values for variables 
70 float mJiorT; // offset from m_h to ht 
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5 float er_th, er_h; // Decision variables 

float ehjh_l, ehjh_r, eh_th_s; // left, right, and right's sign 

float drdth_llnc; // finite difference increments 

float drdhjlnc, drdh_ilnclnc; 
10 float dhdthjlnc; 

float erjilncj", er_hlnc_jr, er_hlnc_h; 

float r"rl; 
float th = thl; 
15 floath = hl; 

#ifdef INCREMENTAL_TRIG 

float prevCos, prevSin; 
#endif 
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lllllllllllllllllllllllllllll LINE DRAWING 111111111111111111111111111111 
llll I III FIRST VOXEL //// 



25 // First voxel is easy when rounding to integer endpoints, but more 

// tricky when using sub-voxel accurate endpoints. Sometimes the voxel 
// the endpoint is actually in should NOT be drawn, because the ideal 
// line would have just nipped the corner of that voxel. What follows 
// is code to determine to back up appropriately so that the first 

30 // voxel drawn by the standard algorithm is the correct one. 

// All the computation done also gives us what we need to know in order 
// to skip regions of the line drawing code that do not apply. 

ra_r = r; 
35 m_h = h; 

m_th = th; 

// Compute partials (or expression derived from partials and 

// implicit function) at current point to determine dominance. 
40 m_cos = cos(m_th * radsPerStep - tht); 

m_sin = sin(m_th*radsPerStep - tht); 

if (m_cos < 0) rn_cos = 0; // clamp to +- 90 degree range 

drdth_r = m__cos; 

drdmj5 = drdth_r<0 ?-l : 1; 
45 drdthj - m_r * m_sin * radsPerStep * drdth_s; 

drdth_r *= drdth_s; 



dhdth_r = m_cos*m_cos; 
dhdthj = radsPerStep*xi; 

drdhj = one_m_eta * m_r*m_r - rt2; 



// determine dominance rankings 

enum { R_DOM_H = 0x4, R_DOM_TH = 0x2, H_DOM_TH = 0x1 }; 
55 unsigned char bits = 

// ((fabs(drdhj) > fabs(drdhj-) || hlgnore) « 2) | // r dominates h 
((drdhj > 0 1| hlgnore) « 2) | // r dominates h 
((fabs(drdthj) > fabs(drdth_r)) « 1) | II x dominates th 
((fabs(dhdthj) > fabs(dhdth_r)) || thlgnore « 0); // h dominates th 
60 if (thl=th2) bits |= (H_DOM_TH|R_JDOM_TH); // can't be TH-dominant 

/* 

switch (bits) { 

case 0x1 :// 001 r<h r<th h>th => h-th, h-dom 
65 case 0x0: // 000 r<h r<th h<th => h-th, th-dom 
break; 

case 0x3: //Oil r<h r>th h>th => r-h, h-dom 
case 0x7: //ill r>h r>th h>th => r-h, r-dom 
break; 

70 case 0x6: // 1 10 r>h r>th h<th => r-th, r-dom 
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5 case 0x4:// 100 r>h Kth h<th => r-th, th-dom 
break; 

// case 0x2: // 0 10 r<h r>th h<th => inconsistent 
// case 0x5: // 101 r>h r<th h>th => inconsistent 

10 } 
*/ 

if(r<l) 
{ 

// IGNORE WINDUP TO AVOID BACKING UP TO r<0 

15 

// IF ETA BIG, THEN R-DOMINANCE IS IMPOSSIBLE 
if (eta>1.0) bits ~(R_DOM_H); 

// IF ACTUAL START POINT IS VERY CLOSE TO R-0, THEN LINE COULD BE 
20 // TH-DOM EVEN THOUGH IT IS NOT ACCORDING TO SLOPE CHECK AT R=0.5. 

// JUMPING DIRECTLY TO REGIONH_PLUS IS WRONG IN THAT CASE, SO ACT AS IF 
// THETA IS DOMINANT. IF IT REALLY ISNT TH-DOM, THEN WE'VE ONLY DONE 
// A BIT OF EXTRA WORK. 

if ((bits & HJX)M_TH) && rlf<0.5) bits ~(H_DOM JTH) ; 

25 

setVoxeI(r,th,h,COLOR_CYAN); 

} 

else //(r>=l) 
{ 

30 // BASED ON DOMINANCE, BACK UP TO REACH A POINT 

// AT WHICH DRAWING OF THE FIRST VOXEL CAN BE DECIDED 
// LIKE ANY OTHER MID-LINE VOXEL DECISION 
if ((bits & (R_D OM_H|R_DOM_TH)) = 0x0) { // h-th 
unsigned char hDom = (bits & H_DOM_TH); // h dominant? 
3 5 bool belowMidpoint = hDom 

? hSgn*zl < hSgn*mJi 
:thlf<m_th; 
if (hDom || belowMidpoint) h -= hSgn; 
if(m_r»=0){ 
40 if (!hDom || ! belowMidpoint) th -= 1; 

} 

} 

else if ((bits & (R_DOM_TH|H_DOM TH)) = (R DOM_TH|HJDOM_TH)) { // r-h 
if (eta>l .0) bits &= «-(R_DOM_H); 
45 unsigned char rDom = (bits & R_DOM_H); // r dominant? 

int sgn = rDom 
? start_minus? -1 : 1 
: (hSgn^y^start^minus ? -1 : 1; 
bool belowMidpoint = rDom 
50 ? sgn*rlf < sgn*m_r 

: sgn*zl < sgn*m_h; ^ 
if (rDom I| (beIowMidpoint&&polarL2J= : 0)) r -= start_minus?- 1:1; r, - — — ~ L ' 
if (IrDom || IbelowMidpoint) h -= hSgn; 

5 5 else if ((bits & (R_DOM_H|H_DOMjra)) == R_DOM_H 

|| (hlgnore&&thlgnore)) { // r-th 

unsigned char rDom = (bits & R_DOM_TH); // r dominant? 

int sgn = start_minus ? -1 : 1; 

bool belowMidpoint » (sgn *r If < sgn*m_r); 
60 if (rDom || belowMidpoint) r -= sgn; 

if (IrDom || IbelowMidpoint) th — 1; 

} 

if (m_th!=th) // we stepped back in theta, recompute sin & cos 
65 { 

m_th = th; 
#ifdef INCREMENTAL_TRIG 

prevCos = m_cos; prevSin = m_sin; 
// Backwards step — note change of signs for bW! 
70 m_cos =* prevCos - (aW*prevCos - b W*prevSin); 
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m_sin = prevSin - (aW*prevSin + bW*prevCos); 
#eise 

m_cos = cos(m_th*radsPerStep - tht); 
m_sin =* sin(m_th*radsPerStep - tht); 
#endif 
} 

}//endif(r>-l) 

// ANY OTHER CASE MUST BE BORDERLINE CASES THAT ARE VERY CLOSE TO BEING 
// DOMINANT IN TWO OR MORE VARIABLES. IN THESE CASES, DO NOT BACKUP. 



// DISTRIBUTOR CAP: JUMP TO APPROPRIATE STARTING REGION 

bits |= ((start_minus&&r I !=0) « 3); 

switch(bits) 

{ 

// MUST DO ALL INITIALIZATION EXPECTED BY TARGET REGION BEFORE JUMPING 

case 016: // 1 1 10 r>h r>th h<th => r-th, r-dom 
case 017: // 1 1 1 1 r>h r>th h>th => r-h, r-dom 

goto REGIONJfcJVONUS; 

break; 

case 01 1 : // 1001 Kh r<th h>th => h-th, h-dom 
case 013: // 101 1 r<h r>th h>th «> r-h, h-dom 
m_th = th + 0.5; 

TRIG_HALF_STEP(m_cos,m_sin); 
goto R£GION_H_MINUS; 
break; 

case 010: // 1000 r<h r<th h<th => h-th, th-dom 
case 014: // 1 100 r>h Kth h<th => r-th, th-dom 
m_th = th+ 1; 

TRIG_WHOLE_STEP(m_cos,m_sin); 
goto REGION_THETA_MINUS ; 
break; 

case 000: // 0000 r<h Kth h<th => h-th, th-dom 
case 004: // 0100 r>h Kth h<th => r-th, th-dom 
m_th = th+ 1; 

TRIG_WHOLE_STEP(m_cos,m_sin); 
goto REGION_THETA_PLUS; 
break; 

case 001:// 0001 Kh Kth h>m => h-th, h-dom 
case 003:// 0011 Kh r>th h>th => r-h, h-dom 
m_th =■ th + 0.5; 

TRIG__HALF_STEP(m_cos,m_sin); 

goto REGIONHJPLUS; 

break; 

case 006: // 01 10 r>h r>th h<th => r-th, r-dom 
case 007:// 01 11 r>h r>th h>th r-h, r-dom 
m_th = th + 0.5; 

TRIG_HALF_STEP(m_cos,m_sin); 

goto REGION_R_PLUS; 

break; 

} 



//////// REGION R- ////// 

REGION_R_MINUS: 

m_r = r-I; 

m_th = th+0.5; 

m_hoff= h+hHalflnc-ht; 

TRIG_HALF_STEP(m_cos,m_sin); 

// The solution in h for the simultaneous equations: 

// r-h slope = 1 && F(r,h) = 0. 

// More robust than just r-h slope test when r-h hyperbola's 

// asymptotes have slope very near +-1 or +-infinity 
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5 drdhj = one_rn_eta*rnj-*m_r - rt2; 
drdhjlnc = one_ni_eta*( 1 -2*m_r); 
drdhjlnclnc = 2*one_m_eta; 

drdthjr = mjcos; 
10 drdthjs - drdth_r < 0 ? 0 : 1; 

drdthjlnc = m_sin * radsPerStep * drdth_s; 
drdthj = m_r * drdthjinc; 
drdth_r *= drdthjs; 

1 5 erjh = m_r * m_cos - rt; 

erji = eta*(m_r*m_r - rt2) - m_hoff*rnJioff, 
er_hlnc__r « eta*(l-2*mjr); 
erjilncjr = 2*eta; 
20 er_hlnc_h = -hSgn*2*m_ho£f - 1; 
// erjilncjih = -2; 

while((hIgnore || drdhj > 0) && 

(thlgnore || drdth_s=0 || drdthj < -drdthj ) && 
25 r>rmin 
) 

{ 

m r-= 1; 

30 

if (th<th2 && erjh < 0) { 
th+=M; 
mJh+= 1; 

TRIG WHOLE STEP(m cos,m sin); 

35 

drdthj* = m_cos; 
drdthj; = drdth_r < 0 ? 0 : 1; 
drdthjlnc = m_sin * radsPerStep * drdth_s; 
drdthj = m_r * drdth_llnc; 
40 drdthj- *= drdthjs; 

erjh = m_r * m cos - rt; 

} 

else { 
er th -= m cos; 
45 } " 

if (h*hSgn<hmax && erji < 0) { 
h += hSgn; 
drdthj — drdthjlnc; 
50 erji += erjilncji; 

erjhlncjh -= 2; // erjilncjih 

} 

drdhj += drdhjlnc; 
5 5 drdhjlnc += drdh JIncInc; 

erji += er_hlnc_r; 
er__hlnc_r += erjilnc_rr; 

60 setVoxel(r,th,h,COLOR_YELLOW); 

} 

if (r<=rmin && rmin=r2 && th=th2 && h=h2) 
return; // we're done here 

65 

//////// REGION H- ////// 
REGION_H_MINUS : 

if (Ihlgnore) { 
70 m_r = r-0.5; 
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5 mjioff =* h+hSgn-ht; 

drdhj = mjioff; // useful in nearly vertical lines to 

// make sure we don't cross into the PLUS regions 

1 0 dhdthj = m JiofPrnJioff; 

dhdthJInc = hSgn*2*m_hoff+ 1; 
dhdth_r = xi * (stepsPerRad*hSgn-xi); 

er_h » eta*(m_r*m_r - rt2) - dhdthj; 
1 5 er_hlnc_r = eta*(l -2*m_r); 
erJiIncjT = 2* eta; 

eh_thj = xi * m_sin; 
eh_th_r = m_cos; 
20 ehjh_s « hSgn*(ehJhj- < 0 ? -1 : 1); 
ehjh_r *= mjioff; 

while(hSgn*drdhJ<=0 && 
(thlgnore || 
25 dhdthj>dhdth_r 
//hSgn*h<dhdth_r 
|| xi=0)&& 
hSgn*h<hmax) 

{ 

30 h+=hSgn; 

mjioff -h= hSgn; 

if (th<th2 && ehjh_s*ehjhj < ehjh_j>*eh th r) { 
th+= 1; 
35 mjh+=l; 

TRIG_WHOLE_STEP(m_cos,m__sin); 

eh_th_l = xi * m sin; 
eh Jh_r = m_cos; 
40 eh Jh_s = hSgn*(eh Jh_r < 0 ? - 1 : 1 ); 
eh_th_r *= mjioff; 
} else { 

eh__th_r += m_cos*hSgn; 



if (erji > 0 && r>rmin) { 
r-=l; 
m_r -= 1; 

er Ji +■ er hlnc_r, 
50 erjilncjr += er_hlnc_rr; 

} 

dhdthj -h= dhdth JInc; 
erh -= dhdth_llnc; 
55 dhdthjlnc += 2; 

drdhj = m_hofT; 

setVoxel(r,th,h,COLORJ3REEN); 

60 } 

} // end if(!hIgnore) 

//////// REGION THETA- /// 

65 

if (Ithlgnore || r<l) { // Can always be theta dominant at polar origin 
mJh+=0.5; 

TRIG_HALF_STEP(m_cos,m_sin); 
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REGION_THETA_MINUS: 

m_r = r - 0.5; 

m_hoff = h+hHalflnc-ht; 

drdth_r = mjcos; 

drdthjs — drdth jr < 0 ? 0 ; 1; // should never be negative 
drdth_J = m_sin*drdth_s; 

while(drdthj 0 && th<th2) 
{ 

eh_thj = xi * m_sin; 
ehjthjr = mjcos; 

eh_th_s = hSgn*(eh_th_r < 0 ? -1 : 1); 
eh_th_r *= mjioff; 

erjh - m_r * m_cos - rt; 

if (hSgn*h<h2*hSgn && ehjh_s*eh_thj > ehjh_s*ehjhj-) { 
h += hSgn; 
m_hofT+= hSgn; 

} 

if (er_th > 0 && r>rmin) { 
m_r — 1; 

} 

th-H= 1; 
m_th+= 1; 

TRIG_WHOLE_STEP(m_cos,m_sin); 
drdth_r = m_cos; 

drdth_s = drdthjr < 0 ? 0 : 1 ; // should never be negative 
drdthj = m_sin*drdth_s; 

setVoxel(r,th,h,COLOR__MAGENTA) ; 

} 

//////// REGION THETA+ /// 
REGION_THETA_PLUS : 

m_r = r+0.5; 

m_hoff= h + hHajflnc - ht; 

dhdthj » radsPerStep*xi; 
dhdth_r = m_cos*m_cos; 

// Implicit fh of just theta. Derived from: 
// dr/dth=I and F(r,th) = 0 

// More robust than plain dr/dth check for lines close to polar origin, 
drdthj - dhdth_r - rt*m_sin*radsPerStep; 

while(drdthj>0 && 

hSgn*dhdthJ < dhdthj- && 
th<th2 && 
m_cos > 0) 

{ 

eh_thj = xi * m_sin; 
eh_thjr = m_cos; 

eh_th_s = hSgn*(eh_th_r < 0 ? -1 : 1); 
ehjthjr *= mjioff; 

er_th = m_r * m_cos - rt; 

if (hSgn*h<h2*hSgn && 

eh_th_s*ehjhj > ehjh_s*eh_th_r) { 
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h += hSgn; 
m_hoff hSgn; 

} 

if(r<r2&& er_th<0) { 
r+= I; 
m_r+= 1; 

} 

th-H= 1; 
m_th+= 1; 

TRIG_WHOLE_STEP(m_cos,m_sin); 

dhdth_r = m_cos*m_cos; 

drdthj = dhdthjr - rt*m__sin*radsPerStep; 

setVoxel(r,th,h,COLOR_RED); 

} 

} // end f thlgnore 

//////// REGION H+ ////// 

m_th - th + 0.5; 
#ifdef INCREMENTAL^TRIG 

// this is a haJf step from the PREVIOUS cos and sin values 
m_cos = prevCos - (aH*prevCos + bH*prevSin); 
m_sin = prevSin - (aH*prevSin - bH*prevCos); 
#else 

m_cos = cos( m_th* radsPerStep - tht ); 
m_sin = sin( m_th*radsPerStep - tht ); 
#endif 

REGION_H_PLUS: 

m_r B r + 0.5; 
mjioff = h+hSgn - ht; 

// The solution in h for the simultaneous equations: 

// r-h slope = 1 && F(r,h) = 0. 

// More robust than just r-h slope test when r-h hyperbola's 

// asymptotes have slope very near +/- 1 or +/-infinity 

drdhj =* one_m_eta*m__r*m_r - rt2; 

drdhjlnc = one_m_eta*( 1 -2*nw); 

drdhjlnclnc = 2*onejn_eta; 

ehjtti J 5=5 xi * mjsin; 
eh_th_r = m_cos; 

ehjh_s - hSgn*(eh_thr < 0 ? -1 : 1); 
eh_th_r *= mjiof?, 

er_h = eta*(m_r*m_r - rt2) - m_hofi*m_hofT; 

erjimc_r== eta*(l+2*m_r); 

erjilncjr = 2*eta; 

erjilncji = -hSgn*2*mJioff - 1; 

while(drdhj<0 ||eta> 1.0) 
{ 

if (h*hSgn >= h2*hSgn) return; // gone far enough-skip region r+; 

h +- hSgn; 
m_hoff-H=hSgn; 

if (th<th2 && ehjth_s*ehjhj < eh th_s*eh_th r) { 
th+»l; 
mJh+= I; 

TRIG_WHOLE_STEP(m_cos, m_sin); 
eh_th_i = xi * m_sin; 
eh th_r = m_cos; 
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eh_th_s - hSgn*(eh_th_r < 0 ? -1 : 1); 
eh_th_r *= m_hoff; 

} 

else { 

eh_th_r += m_cos*hSgn; 

} 

if(r<r2&&erji<0) { 
r+=l; 
m_r-H= 1; 

drdhj += drdh JInc; 
drdh JInc += drdh_llnclnc; 
er_h += er_hlnc_r; 
er_hlnc_j += er__hlnc_rr, 

} 

er_h += er_hlnc_h; 
er_hlnc_h -= 2; 

setVoxel(r,tb,h,COLORj3REEN); 

} 



//////// REGION R+ ////// 
REGION_R_PLUS: 
m_r - r+1; 

m_hofF= h+hHaiflnc-ht; 

er_th = m_r * m_cos - rt; 

er_h » eta*(m_r*m_r - rt2) - m_hoff*m_hoff; 

er_hlnc_r « eta*(l+2*m_r); 

er_hInc_.iT = 2*eta; 

erjilncji = -hSgn*2*m_hoff - 1; 

while(r<r2 && h*hSgn <= h2*hSgn && th<=th2) 
{ 

r+=l; 
m_r-H= 1; 

if (th<th2 && er_th > 0) { 

th+=l; 

m_th+= 1; 
tfifdef INCREMENTAL TRIG 

TRIG_WHOLE_STEP(m_cos,m_sin); 
#else 

m_cos — cos( m_th * radsPerStep - tht ); 
#endif 

er_th = m_r * m_cos - rt; 

} 

else { 
er_th 4— m_cos; 

} 

if (hSgn*h < hSgn*h2 && erji > 0) { 
h += hSgn; 
erjh += er_hlnc_h; 
erjilncji -= 2; 

} 

er_h +=* er_hlnc_r, 
erjhlncjr +— erjilnc_rr; 

setVoxel(r,th,h,COLOR_YELLOW); 

} 

} 
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Having described the invention and a prederred embodiment thereof, we claim: 

CLAIMS 

A method for rendering, on a volumetric display having a plurality of voxels 
defined by cylindrical coordinates, a rasterized line that approximates a desired 
line, said method comprising: 

rendering, on said volumetric display, a first voxel from said plurality of 

voxels, the first voxel being representative of a first point from said desired 
line, said first point being selected from a first constituent line-segment of 
said desired line, said first constituent line-segment having a first dominant 
direction; 

selecting a second point from said desired line; 

determining whether said second point is on a second constituent line-segment 
of said desired line, said second constituent line-segment having a second 
dominant direction different from said first dominant direction; 

on the basis of said determination, rendering a second voxel from said plurality 
of voxels, said second voxel being representative of said second point on 
said desired lineand having a location that minimizes a distance between 
said second voxel and said second point on said desired line. 

The method of claim 1 wherein determining whether said second point is on a 
second constituent line-segment of said desired line comprises 

evaluating a difference between a threshold and a transition-indicating 

expression evaluated at said second point, said difference being indicative 
of a transition from said first constituent line-segment to said second 
constituent line-segment. 

The method of claim 2, further comprising determining said second dominant 
direction. 
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4. The method of claim 3, wherein determining said second dominant direction 
comprises selecting said second dominant direction from a subsequent direction 
in a pre-defined progression of dominant directions. 

5. The method of claim 4, further comprising defining said progression of 
dominant directions to be a sequence consisting of: 

a radial direction followed by 

an axial direction followed by 

an angular direction, followed by 

said axial direction, followed by 

said radial direction. 

6. The method of claim 2, further comprising selecting said transition indicating 
expression to include an expression indicative of a slope of said desired line at 
said second point. 

7. The method of claim 6, wherein said first dominant direction defines first and 
second non-dominant directions and said slope is defined to be a ratio of how 
much said desired line changes in said first dominant direction with respect to a 
change in one of said first and second non-dominant directions. 

8. The method of claim 6, wherein said transition indicating expression is selected 
to identify a line having a point at which said slope is equal to said threshold. 

9. . The method of claim 1, wherein rendering said second voxel comprises 

applying a midpoint method to select said second voxel from a plurality of 
candidate voxels. 

10. A method of rendering, on a display having a plurality of voxels defined by a 
cylindrical coordinate system, a rasterized line that approximates a desired line, 
the method comprising: 
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rendering a first voxel on said rasterized line, said first voxel being 
representative of a first point on said desired line; 

determining, on the basis of a local slope of said desired line and a sign of a 
discriminant of said desired line in a neighborhood of said first point, a 
direction relative to said first voxel of a second voxel on said rasterized 
line; 

rendering said second voxel on said rasterized line, said second voxel being 
displaced from said first voxel in said direction. 

11. The method of claim 10 further comprising evaluating said discriminant at a 
midpoint between two candidate voxels from which said subsequent voxel is to 
be selected. 

12. The method of claim 10 wherein determining a direction comprises partitioning 
said desired line into a plurality of constituent line-segments, each constituent 
line-segment being associated with a range of local slopes of said desired line. 

13. The method of claim 12 further comprising updating said discriminant. 

14. The method of claim 12 wherein updating said discriminant comprises: 

incrementing said discriminant by a selected amount when said subsequent 
voxel is displaced from said first voxel in a first direction; and 

incorporating information indicative of a displacement of said second voxel in 
a second direction. 

15. A volumetric display comprising: 

a display surface rotatable about an axis; 

an voxel-rendering element configured to render a selected voxel on said 
display surface; and 

a rasterizer in communication with said voxel-rendering element, said 

rasterizer being configured to direct said voxel-rendering elemen to render, 
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on said display surface, a rasterized line that approximates a desired line 
by: 

rendering, on said display surface, a first voxel defined by first 
cylindrical coordinates, said first voxel being representative of a 
first point from said desired line, said first point being selected 
from a first constituent line-segment of said desired line, said first 
constituent line-segment having a first dominant direction; 

selecting a second point from said desired line; 

determining whether said second point is on a second constituent line- 
segment of said desired line, said second constituent line-segment 
having a second dominant direction different from said first 
dominant direction; 

on the basis of said determination, rendering a second voxel 

representative of said second point on said desired line, said second 
voxel having a location defined by second cylindrical coordinates 
that minimizes a distance between said second voxel and said 
desired line. 

A method for rendering, on a volumetric display having a plurality of voxels 
defined by cylindrical coordinates, a rasterized line that approximates a desired 
line, said method comprising: 

determining whether first and second points on the desired line are on first and 
second constituent line-segments of said desired line, said first and second 
constituent line-segments having different dominant directions; 

on the basis of said determination, incrementally proceeding from a first voxel 
representative of said first point to a second second voxel representative of 
said second point on said desired line, 

rendering said second voxel at a location defined by cylindrical coordinates 
that minimizes a distance between said second voxel and said desired line. 
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